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Abstract 

We have performed a high-precision Monte Carlo study of the dynamic 
critical behavior of the Swendsen-Wang algorithm for the two-dimensional 3- 
state Potts model. We find that the Li-Sokal bound (r intj £ > const x Ch) is 
almost but not quite sharp. The ratio T mtt s/CH seems to diverge either as a 
small power (~ 0.08) or as a logarithm. 
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1 Introduction 



Monte Carlo (MC) simulations U §, |, g have become a standard and powerful 
tool for gaining new insights into statistical-mechanical systems and lattice field the- 
ories. However, their practical success is severely limited by critical slowing-down: 
the autocorrelation time r — that is, roughly speaking, the time needed to produce 
one "statistically independent" configuration — diverges near a critical point. More 
precisely, for a finite system of linear size L at criticality, we expect a behavior r ~ L z 
for large L. The power z is a dynamic critical exponent, and it depends on both the 
system and the algorithm. 

Single-site MC algorithms (such as single-site Metropolis or heat bath) have a 
dynamic critical exponent z > 2. This makes it very hard to get high-precision data 
very close to the critical point on large lattices. 

In some cases, a much better dynamical behavior is obtained by allowing non-local 
moves, such as cluster flips.[] The Swendsen-Wang (SW) cluster algorithm [^] for the 
g-state ferromagnetic Potts model achieves a significant reduction in z compared to 
the local algorithms: one has z between and ~ 1, where the exact value depends 
on q and on the dimensionality of the lattice ||. The most favorable case is the 
two-dimensional (2D) Ising model (q = 2), in which z < 0.3 |7|, [|, |], [l(| (see below). 
In other cases, the performance of the SW algorithm is less impressive (though still 



quite good): e.g., z = 0.55 ± 0.03 for the 2D 3-state Potts model fill , z ~ 1 for the 
2D 4-state Potts model JTTJ 0, |T3} [TJ, and z « 1 for the 4D Ising model |I§ |TB . 
Clearly, we would like to understand why this algorithm works so well in some cases 
and not in others. We hope in this way to obtain new insights into the dynamics of 
non-local Monte Carlo algorithms, with the ultimate aim of devising new and more 
efficient algorithms. 

There is at present no adequate theory for predicting the dynamic critical behavior 
of an SW-type algorithm. However, there is one rigorous lower bound on z due to 



Li and Sokal [IT]. The autocorrelation times of the standard (multi-cluster) SW 



algorithm for the ferromagnetic g-state Potts model are bounded below by a multiple 
of the specific heat: 

Tmt,M, 7int,£, 7"cxp > Const X C H ■ (1-1) 

Here N is the bond density in the SW algorithm, £ is the energy, and Ch is the 
specific heat; r int and r exp denote the integrated and exponential autocorrelation 
times, respectively [|, §]. As a result one has 

a 

where a and v are the standard static critical exponents. Thus, the SW algorithm for 
the g-state Potts model cannot completely eliminate the critical slowing-down in any 



model in which the specific heat is divergent at criticality. The bound (|1 . 1|) / (|1 . 2|) has 
also been proven to hold |13[] for the direct SW-type algorithm [I7| for the Ashkin- 
Teller (AT) model p], HI- 



1 See H, [|, H for reviews of collective-mode Monte Carlo methods. 



2 



The important question is the following: Is the Li-Sokal bound (1.1)/ ( |OD sharp 



or not? An affirmative answer would imply that we could use the bound to predict the 
dynamic critical exponent (s) z given only the static critical exponents of the system. 
There are three possibilities: 



i) The bound ( |1.1| ) is sharp (i.e., the ratio t/Ch is bounded), so that ( |1.2| ) is an 
equality. 

ii) The bound is sharp modulo a logarithm (i.e., t/Ch ~ log p L for p > 0). 

hi) The bound is not sharp (i.e., t/Ch ~ LP for p > 0), so that ( |L2| ) is a strict 
inequality. 

Unfortunately, the empirical situation, even for the simplest cases, is far from 
clear. Let us review the status of this problem for the 2D Potts models. For the 
Ising case, the numerical data J71 |8], |], |It| are consistent with a power-law behavior 
with z- int£ ranging from 0.35 ± 0.01 to 0.25 ± 0.01 || |10[| . However, the data are 



also consistent with a logarithmic behavior z- mt ,e = x log ||. In [|T3[] we reanalyzed 
the high-precision data of Baillie and Coddington 0, [1O|0, and we found that the 



ratio T mti e/CH behaves most likely as a pure power law ~ LP with p = 0.060 ± 0.004 
(statistical error only) or as a logarithm ~ A + BlogL. This means that the bound 
( |1.1|) is either non-sharp by a small power, or else it is sharp modulo a logarithm (in 
the latter case, the leading term would be r int: £ ~ log 2 L). It is extremely difficult to 
distinguish between these two scenarios with lattice sizes up to only L = 512. 

The 3-state Potts model was first considered under this perspective in [II]]. The 
dynamic critical exponent was found to be Zint,£ — 0.55 ± 0.03, which is significantly 
larger than the exact result a/v = | = 0.4 [2~I]. So, it seemed that the bound ( |1.2| ) 
was not sharp at all. 

The 4-state Potts model is rather peculiar: the naive fit to the data, z- mt ,e = 0.89 ± 
0.05 PHI , is smaller than the (exactly known) value of a/v = 1 [22]. The explanation 



of this paradox is that the true leading term in the specific heat has a multiplicative 
logarithmic correction, Ch ~ Llog~ 3 / 2 L |14|, |23|, B3, 5HJ. Indeed, a naive power- law 



fit to the specific heat yielded a/v = 0.75 ± 0.01, consistent with the bound (|1.2|) . 
A high-precision study of this model was carried out in [ 12] , |T3|.p| Naive power- law 



fits to the data showed that the bound ( |1.2|) was satisfied: z intt £ = 0.876 ± 0.012 
and a/v = 0.768 ± 0.009. On the other hand, the behavior of the ratio t^s/Ch 
was consistent only with two scenarios: a power law ~ LP with p = 0.118 ± 0.012 



(statistical error only), or a logarithmic growth ~ A + B log L. This means that (|Q 



fails to be sharp either by a small power or by a logarithm. In conclusion, there are 



2 The exact value of the specific heat at finite L was taken from the paper by Ferdinand and 
Fisher @. 

3 In jl2l the "embedding" version of an SW-type algorithm for the AT model was used. This 
algorithm reduces to the standard SW algorithm at the Ising subspace, but not at the 4-state Potts 
subspace. However, it was shown numerically that both algorithms belong (as expected) to the same 
dynamic universality class at the 4-state Potts subspace. 
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two likely behaviors for the autocorrelation time: either 7i n t,£ ~ L q log 3 ^ 2 L with 

-| in 

q w 1.12, or else r- mt ^ ~ I/log ' L. In both cases, we find multiplicative logarithmic 
corrections to the autocorrelation time, which make the numerical analysis extremely 
difficult. These two scenarios for the ratio t^^/Ch coincide with those obtained for 
the Ising case. 

In summary, for the 2D Potts models with q = 2 and q — 4, the Li-Sokal bound 
might be either sharp modulo a logarithm or else not sharp with a very small power 
p = z — ajv (0.04 < p < 0.12). Moreover, if we interpolate between these models 
by following the self-dual (critical) curve of the AT model |Tj| |T9|] , we obtain [T2|, [T3 



the same two scenarios: the ratio t^^/Ch either grows like a power law with a very 
small power p (which increases slightly as we move from the Ising model to the 4- 
state Potts model along the AT self-dual curve), or else grows like a logarithm. Thus, 
there is some kind of continuity along the AT self-dual curve for the dynamic critical 
behavior of the SW algorithm.^ 

There results have led us to reappraise the status of the Li-Sokal bound for the 
2D 3-state Potts model. In |TT[ the bound was declared not sharp on the basis 
of numerical evidence suggesting that p = z — a/u = 0.15 ± 0.03. However, this 
value is not much larger than that obtained [12, 13] for the 4-state Potts model; this 
suggests that the data for the 3-state model might also be consistent with a logarithm. 
Indeed, a closer look at the results of [|ll|] concerning the 3-state model reveals that 
the lattices studied were not very large (L < 256), and the statistics were not very 
high (the number of iterations for the L = 256 lattice was at most 25 x 10 3 T int £). 
This motivated us to reconsider this model and extend the results of |TT| to larger 
lattices, with higher statistics. 

This paper is organized as follows: Section |2] reviews the basics of the Swendsen- 
Wang algorithm for the Potts models, as well as the proof of the Li-Sokal bound for 
these models. In Section |3| we describe our Monte Carlo simulations. In Section ^ 
we discuss in detail our methods of statistical data analysis. Finally, in Section |5] 
we present the analysis of our numerical (static and dynamic) data, culminating in a 
discussion of the sharpness of the Li-Sokal bound. 



2 Basic set-up and notation 

2.1 Potts model and Swendsen— Wang algorithm 

The g-state Potts model assigns to each lattice site i a spin variable Oi taking 
values in the set {1, 2, ... , q\\ these spins interact through the reduced Hamiltonian 

ftpotts = -/^EOWi- 1 )' (2- 1 ) 

(ij) 

where the sum runs over all the nearest-neighbor pairs (ij). To simplify the notation 
we shall henceforth write S aulJ] = 5 ab for a bond b = (ij). The ferromagnetic case 

4 A similar study was carried out in [[r7| for the single-cluster version of the algorithm, but the 
lattices were not very large (L < 256). 



4 



corresponds to (3 > 0. The partition function is defined as 



£exp 

M 



/3£(<^-i) 



Finally, the Boltzmann weight of a configuration {a} is given by 



WWW) 



n 



,P(Sa b -l) 



H(l-p + P 5 ab ) 



(2.2) 



(2.3) 



where p = 1 — e _/3 . 

The idea behind the Swendsen-Wang algorithm 0, |6], is to decompose the 
Boltzmann weight by introducing new dynamical variables n& = 0, 1 (living on the 
bonds of the lattice), and to simulate the joint model of old and new variables by 
alternately updating one set of variables conditional on the other set. The Boltzmann 
weight of the joint model is 



^•oint({0-},{n}) = ~ p)S nb ,0 + pSa b S nb ,l\ 

L b 



(2.4) 



The marginal distribution of ( |2.4| ) with respect to the spin variables reproduces the 
Potts-model Boltzmann weight ( |2.3|) . The marginal distribution of ( |2.4|) with respect 
to the bond variables is the Fortuin-Kasteleyn [p7| , |2"5 , |29| random-cluster model with 
parameter q: 



W RC ({n}) = \ 



C({n}) 



(2.5) 



n p n v 

b: n b =l b: nt,=0 

where C({n}) is the number of connected components (including one-site components) 
in the graph whose edges are the bonds with n& = 1. 

We can also consider the conditional probabilities of the joint distribution (|2.4j) . 
The conditional distribution of the {n} given the {a} is as follows: Independently 
for each bond b = (ij), one sets = when Oi ^ <jj, and sets = and 1 with 
probabilities 1 — p and p when Oi = <jj. Finally, the conditional distribution of the 
{a} given the {n} is as follows: Independently for each connected cluster, one sets all 
the spins Oi in that cluster equal to the same value, chosen with uniform probability 
from the set {1, 2, ... , q}. 

The Swendsen-Wang algorithm simulates the joint probability distribution (|2.4| ) 
by alternately applying the two conditional distributions just described. That is, we 
first erase the current {n} configuration, and generate a new {n} configuration from 
the conditional distribution given {a}; we then erase the current {a} configuration, 
and generate a new {a} configuration from the conditional distribution given {n}. 



2.2 Li-Sokal bound 

To prove the Li-Sokal bound we first notice that the transition matrix Psw of the 
Swendsen-Wang algorithm can be written as a product 



sw 



-Pbond -Pspin 



(2.6) 
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where Pbond (the update of the bond variables) and P sp in (the update of the spin 
variables) are given by the conditional expectation operators E( ■ |{cr}) and E( ■ \ {n}), 
respectively. 

The strategy behind the proof is to compute explicitly the autocorrelation function 
at time lags and 1 for a suitable observable O (which we will choose to be a "slow 
mode" of the algorithm). The unnormalized autocorrelation function is defined as 

Coo{t) = (0(s)0(s + t))-(0) 2 (2.7) 

(where the expectations are taken in equilibrium), and the normalized autocorrelation 
function as 

Pooit) s ^ = ^ . (2.8) 
Coo(0) var(C) 

Then, using some general properties of reversible Markov chains, we will deduce lower 
bounds for the autocorrelation times Tint,© and T exP) e>. These will in turn imply lower 
bounds on the dynamic critical exponents 2i n t,o an d z exPt o- 
For the observable O, we shall use the bond occupation 

Af = Y, n b ■ (2-9) 

b 



From the joint Boltzmann weight ( |2.4j ) it is easy to compute the following bond 
expectation values conditional on the spin configuration {a}: 

E{n b \{a}) = p5 Ub (2.10a) 

From these equations it is easy to compute the mean values (J\f) and (A/" 2 ), and thus 
Cjw(0) = var(jV) = (A/" 2 ) - (A/") 2 : 

(A/") = p(£) (2.11a) 

(A/" 2 ) = p 2 (£ 2 )+p(l-p)(£) (2.11b) 

CW(0) = var(AT) = p 2 var(£) + p{\ - p){£) (2.11c) 

where the energy is defined as 

£ = T,< ■ ( 2 - 12 ) 

6 

The unnormalized autocorrelation function at time lag 1 is given by 

CW(1) = (W(O)jV(l)) - (AT) 2 = var(P bond AT) = var(P(AT|{a})) . (2.13) 
Now, PbondA/" is equal to 

PbondAT ee E(Af\{a}) = P J2K = pS. (2.14) 

b 
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Therefore, 

CW(1) = P 2 var(£) (2.15) 

and 

WM J - CW(0) " ^ + (1 - p)i5 - P C H + (1 - p)i5 ' {2 - lb) 
where the energy density E and the specific heat Ch have been defined as 

E = i<£) (2.17) 
= ^var(£) ee 1 [<£ 2 > - (S) 2 ] , (2.18) 

and V is the number of lattice sites. 

The correlation functions of Af under Psw = -Pbond-Pspin are the same as under the 
positive-semidefinite self-adjoint operator P' sw = P S pin-fbond-Ps P in- This implies (see 
e.g., ||) that we have a spectral representation 

PNN {t) = f l A 1 ' 1 dv{\) (2.19) 
Jo 

with a positive measure dv. From this spectral representation we conclude that 

PMM{t) > PATAr(l) 1 ' 1 • (2.20) 

If we now recall the definitions of the integrated and exponential autocorrelation 
times 

oo 

T"int,AT = - AAW(^) (2-21) 

1 t=-oo 



\t\-^o log PMM{t) 



7"exp,Af = ,11111 T7Z - — (2.22) 



we conclude from ( |2.16| ) / ( |2.20| ) that 



rmtx > \ x 1 + > const x C H (2.23) 

r eX p,AA > i 7-t > const x C H (2.24) 

These are precisely the bounds ( |1 . 1|) . If we take into account the expected behavior 
close to the critical point of the specific heat and the autocorrelation times, we deduce 
immediately the bounds (|1.2|) . 

Similar bounds hold for the autocorrelation times of the energy £ . This can be 
seen from the fact that PbondA/" = pS, which implies that 

C ee (t) = jT 2 CW(i+l) (2.25) 
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and hence 

Pee® = P^ + V > PNN (t) . (2.26) 
This equation proves that 

7exp,£ = T~exp,Af j (2.27) 

and furthermore allows us to bound Tfot,£ above and below in terms of T- mt) M'- 

, T mt,Af - \ 1 r int,AT , , 

TintX < 7Jnt.fi = — 77^ ~ 7j < ~ 7TT • (2-28) 

If the critical slowing-down is not completely eliminated, we expect the factor pjvov(l) 
to approach 1 (from below) as L — > oo. Moreover, irrespective of the presence or 
absence of critical slowing-down, we expect pmnX) to be bounded away from zero 
as L — » oo. Modulo this very weak hypothesis, Q2.28 ) implies the equality of the 
dynamic critical exponents for the energy and bond-occupation observables: 

z iat,8 — z mt,M ■ (2.29) 

Unfortunately, we do not know how to rigorously rule out "exotic" behaviors, in which 
Pa/0v(1) tends to zero as L — > oo and yet r int ^ diverges because the autocorrelation 
function has an extremely long tail. 

Finally, we can define a new ( "energy-like" ) bond observable: the nearest-neighbor 
connectivity 

£' = £ T &, (2-30) 

b 

where 7& equals 1 if both ends of the bond b belong to the same cluster, and 
otherwise. More generally, the connectivity 7^ can be defined for an arbitrary pair 
i,j of sites: 

J 1 if i is connected to j , > 

lij\\ n i) 1 q if 2 is not connected to j 

The interest in £' stems from the fact that the conditional expectation of 8 ab given 
the bond configuration {n} is essentially 7^: 

E(l - 5 ab \{n}) = i^(l- 7b ). (2.32) 

This implies the following relation between the energy and the connectivity densities: 

where k is the coordination number of the lattice (i.e., the number of nearest neighbors 
of any given site), and the connectivity density E' is defined as 

E> = i<0 . (2.34) 
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Furthermore, ( j2.32| ) tells us that 
P spin (kV/2-£) = P spin ^2(l-S ab ) = i^^(l- 76 ) = ±zl( k Vl2-£') . (2.35) 

h 1 h 1 



This implies that 



and hence 



Cs' £ '(t) = (-^-\ C ££ (t + 1) (2.36) 
P£>e>(t) = M^Lil > p£g{t) . (2.37) 



Thus, the bounds ( p..l|) /(L2|) hold also for the autocorrelation times of £'. Further- 



more, we can obtain bounds analogous to (|2.27|) - (|2.29 ) for the observable £': 



7"exp,£ — T eX p,£' (2.38) 

and 

Tint.g < r int ,£' < mt f^ • (2.39) 
Pse{l) 



Using again Q2.28| ) and ( |2.26| ) we arrive at 



- - ' (2 ' 40) 

If pmm(2) is bounded away from zero as L — ► oo, we can conclude that the dynamic 
critical exponents of A/", £ and £' are all equal: 

^int,AT = Z- mtj e = Zint,£> ■ (2-41) 

3 Description of the simulations 

3.1 Observables to be measured 

We have made simulations of the two-dimensional 3-state Potts model at critical- 
ly, 

p = (3 C = \og(l + V3) « 1.00505254, (3.1) 

on a periodic square lattice of linear size L. 

We have measured five basic observables. Three of them have been already de- 
fined: the energy £ (|2.12|) , the bond occupation M Q2.9| ), and the nearest-neighbor 



connectivity £' ( |2.30| ). The other two are 

M 2 = (5>*1 (3.2a) 



X 



q * f \ 2 V 2 



1 a =l \ x 



(3.2b) 
q - 1 
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and 



2mx\/L 




+ 



E 



cr T e 



2mx2/L 



+ 



(3.3a) 



(3.3b) 



where <r x G 
this means <r x 



is the Potts spin in the hypertetrahedral representation [for q = 3 
(cos(27r<T x /3), sin(27TO" a: /3))], V" = L 2 is the number of lattice sites, 
and (x\,X2) are the Cartesian coordinates of point x. The observable T can be 
regarded as the square of the Fourier transform of the spin variable at the smallest 
allowed non-zero momenta [i.e., (±27r/L, 0) and (0, ±2ir/L) for the square lattice]. It 
is normalized to be comparable to its zero- momentum analogue M 2 . 

From these observables we compute the following expectation values: the energy 
density E (|2.17| ), the specific heat Ch ( |2.18|) , the connectivity density E' (|2.34|) , the 
bond density 



N 



V 



(AO 



the magnetic susceptibility 



X 



> 2 > 

the correlation function at momentum (2tt/L, 0) 

1 



V 

and the second-moment correlation length 



1/2 



(3.4) 



(3.5) 



(3.6) 



(3.7) 



2sin(7r/L) 

This definition of the correlation length is not equal to the exponential correlation 
length (=l/mass gap); but it is expected that both correlation lengths scale in the 
same way as we approach the critical point. 

Remark: As a check we have also computed the mean-square size of the clusters, 



So 



c 



(3.8) 



where the sum is over all the clusters C of activated bonds (i.e., with n& = 1) and 
is the number of sites of the cluster C. Using the Fortuin-Kasteleyn identities 
H [27|> |8|, |9f, i1] is not difficult to show that 



(S 2 ) = (M 2 ) . 



(3.9) 



For each observable O discussed above we have measured its autocorrelation func- 
tion poo{t) fl2.8p ; and from this function we have estimated the corresponding inte- 
grated autocorrelation time r intj o Q2-21J) . In Section [| we explain in detail how we 
derived estimates of the mean values and the error bars for both static and dynamic 
quantities. 
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3.2 Summary of the simulations 



We have run our Monte Carlo program on lattices with L ranging from 4 to 1024. 
In all cases the initial configuration was random, and we have discarded the first 10 5 
iterations to allow the system to reach equilibrium; this discard interval is always 
greater than 10 rwf. Q For 4 < L < 256 the total run length was approximately 
10 6 r intj£ ; for L = 512, it was 2.2 x 10 5 r intj£ ; and for L = 1024, it was 6.7 x 10 4 r intj£ . 
In all cases, the statistics are high enough to permit a high accuracy in our estimates 
of the static (error ~ 0.1-0.5%) and dynamic (error ~ 0.5-2%) quantities. 

For 4 < L < 128 our data were obtained from a single long run. For L = 256 
we made two independent runs (with different random-number-generator seeds); for 
L = 512 we made three independent runs; and for L = 1024 we made four independent 
runs. In each case, we discarded the first 10 5 iterations of each run. The individual 
runs are all of length > 10 4 T int £, which is long enough to allow a good determination 
of the dynamic quantities. 

To test the program, we compared the static results on the 4x4 lattice to the exact 
solution (obtained by enumerating all the possible configurations on this lattice). The 
agreement was excellent (see Table [I]). In addition, for all lattice sizes we checked 
the relations (|2.11a| ), (|2.33| ) and ( |3.9| ) between the mean values of static observables; 
the relations ( |2.25[ )/( p.36| ) among the autocorrelation functions C^v, Cg£ and Csi£>\ 
and the relation (|2.16|) between the autocorrelation function pmm^X) an d the static 
observables E and Ch- In all cases the agreement was also good. 

The CPU time required by our program is approximately 7.2 L 2 /is/iteration on an 
IBM RS-6000/370, and 3.6 L 2 /is/iteration on a single processor of an IBM RS/6000 
SP2. The total CPU time used in the project was approximately 1.5 years on an IBM 
RS/6000 SP2. The smallest lattices were run on a IBM RS-6000/370 at NYU, while 
the largest ones (L > 128) were run on the IBM SP2 cluster at the Cornell Theory 
Center. 



4 Statistical analysis of the Monte Carlo data 

In this paper we are aiming at extremely high precision for both static and dynamic 
quantities; and furthermore we need to disentangle the effects of statistical errors from 
the effects of systematic errors due to corrections to scaling. For this, it is essential 
to obtain accurate estimates not only of the static and dynamic quantities of interest, 
but also of their error bars: in this way we will be able (see Section [5]) to perform \ 2 
tests which provide an objective measure of the goodness of fit in each scaling Ansatz. 

In this section we discuss in some detail how we performed the statistical anal- 
ysis of our raw Monte Carlo data. In particular, we describe how to compute the 
estimators for the mean value and the variance of both static and dynamic quanti- 



ties. These methods are based on well-known results of time-series analysis [|3C], |3T 



5 We expect that Ti nt ,£/T exPi £ rj 0.96 and t CXP! £ = r oxp for this algorithm |13[ (see Section p.5[ ) 
So the discard interval is always greater than 10 3 r cxp , which is more than sufficient. 
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which we review briefly in Section f4.1| .p| Then (Section |4.2|) we describe an alternative 
analysis method, based on independent "bunches" , and report the results of detailed 
cross-checks that confirm (with one slight exception) the validity and reliability of 
the standard time-series-analysis method. 



4.1 "Standard" time-series-analysis method 

Let us consider a generic observable O, whose mean is equal to \iq- Its cor- 
responding unnormalized and normalized autocorrelation functions are denoted by 
C 00 (t) = (O(O)C(t)) - (O) 2 and Poo (t) = C O o(t)/C O o(0), respectively. We also 
define the integrated autocorrelation time 

oo 

Tint,0 = - E POoit) • (4-1) 
z t=-oo 

Given a sequence of n Monte Carlo measurements of the observable O — call 
them {Oi, ... , O n } — the natural estimator of the mean \iq is the sample mean 

_ 1 n 

O = (4.2) 

This estimator is unbiased and has a variance 

_ 1 n 

var(O) = — ]T Cao(r - s) (4.3a) 
E ( l -^) C oo{t) (4.3b) 

1 



t=— (n— 1) 



2r in t,ci Coo(0) for n > T iat ,o (4.3c) 
n 

This means that the variance is a factor 2Ti ntj e> larger than it would be if the mea- 
surements were uncorrelated. It is, therefore, very important to estimate the auto- 
correlation time Ti nt) o in order to ensure a correct determination of the error bar on 
the (static) quantity /ie>. 

The natural estimator for the unnormalized autocorrelation function Coo(t) is 

i n-\t\ 

Coo{t) = E (°i ~ Vo)(O i+lt \ - M (4.4) 

n r\ i=i 



if the mean \iq is known, and 



I n-\t\ 



Coo(t) = E (°i ~ 0)(Oi+\*\ ~ °) (4-5) 

n l r l i=i 



6 A review of time-series-analysis methods as applied to MC simulations can be found in |3S 
Appendix C] and in j|. 
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if the mean //© is unknown. We emphasize that, for each t, the estimators Coo{t) 
and C 00(f) are random variables [in contrast to Coo{t), which is a number]. The 
estimator Coo{t) is unbiased, and Coo{t) is biased by terms of order 1/n. The 
covariance matrices of Coo and Coo are the same to leading order in the large- n 



limit (i.e., n ^> r int C i), and we have |30], |31 



J oo 

cov (Coo (t), Coo (u)) - - V [Coo (m) Coo (m + u - t) + C o(m + «) Coo (m - t) 

n ■ 

l = — OC 

+n(t, m, m + u)\ + o ^— ^ , (4.6) 
where i, w > and k is the connected 4-point autocorrelation function 

n(r,s,t) = ({Oi - }io){Oi+r - no)(O i+s - fi ){Oi+ t - fio)) 
-C o(r)Coo(t - s) - Coo(s)C o(t - r) 
-C o(t)Coo(s-r) . (4.7) 

The natural estimator for the normalized autocorrelation function poo{t) is 

Poo(t) = -p— - (4.8) 



if the mean po is known, and 

Poo(*) = ~ ( 4 - 9 ) 

Coo(0) 

if the mean \xq is unknown. The estimators pooif) and Poo{t) are biased by terms of 
order 1/n, as a result of the ratios of random variables in ( [4.81) / Q4.9|) . The covariance 
matrices of poo and p 00 are the same to leading order in 1/n. If the process is 
Gaussian, this covariance matrix is given in the large-n limit by |3lf 



oo 

cov(poo(t), poo{u)) = - [poo(m)poo('rn + t-u) + poo{m + u)p o{m - t) 

^ m=— oo 

+2poo(t)poo(u)p 2 00 (m) - 2poo(t)poo(m)poo(m - u) 
-2poo(u)poo{m)poo{m-t))+o(^j . (4.10) 

If the process is not Gaussian, then there are additional terms proportional to the 
fourth cumulant n(m,t,t — u); these terms are, like those in Q4.10D , of order 1/n. 
The simplest assumption is to consider the stochastic process to be "not too far from 
Gaussian", and drop all the terms involving k. If this assumption is not justified, 
then we are introducing a bias in the estimate of this covariance. 

Finally, we shall take the estimator for the integrated autocorrelation time to be 

1 m 

Tint,© = - E POO® (4-11) 
/ t=-M 
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[or the same thing with poo{t)} where M is a suitably chosen number. The reason 
behind the cutoff M is the following: if we were to make the "obvious" choice M = 
n + 1, then the resulting estimator would have a variance of order 1 even in the limit 
n — > oo; this is because the terms poo{t) with large t have errors (of order 1/n) that 
do not vanish as t grows [cf. ( |4.10| )1, and their number is also large (~ n). Taking 
M<n restores the good behavior of the estimator as n — > oo. The bias introduced 
by this rectangular cutoffQ is given by 

bias(f intj0 ) = ~ Poo(t)+o(-) . (4.12) 

1 \t\>M Vn/ 



The variance of the estimator f int q can be computed from the covariance ( 4.10| ); the 



final result is [32 



_ . 2(2M + 1) 2 
var(r intj0 ) « r int , (4.13) 

where the approximation T- m t,o <M<n has been made. A good (self-consistent) 
choice of M is the following ]32|]: let M be the smallest integer such that M > 



CTi n t : o(M) , where c is a suitable constant. If the normalized autocorrelation function 
is roughly a pure exponential^, then the choice c ~ 6 is reasonable. Indeed, if we take 
Poo(t) = e _<//r and minimize the mean-square error 

MSE(f inti0 ) = bias(f inti0 ) 2 + var(f intj o) (4.14) 

using (|4.12|) / (|4.13|) , we find that the optimal window width is 



M o P t = 2 l0g (27 J ~ 1 ■ (415) 

For njr ~ 10 6 (resp. 10 4 ), we have M opt /r w 6.56 (resp. 4.26). 

As noted above, we expect the estimator f int C i to have a bias of order r intt0 /n, due 
to the nonlinearities in ( |4.8D / (|4.9D .P| To make this bias negligible we need long runs. 
It has been shown empirically that this procedure works fairly well when n > 10 4 ?i nt o 
I- 

Remarks: 1. The estimation of the error bar for the specific heat is a little bit more 
complicated. One can obtain var(C#) by computing var(£), var(£ 2 ), and cov(£,£ 2 ). 
This procedure has a numerical drawback: sometimes the covariance matrix for the 
observables 8 and S 2 is nearly singular; then a small statistical fluctuation can cause 
the estimator of this matrix to be non-positive-definite. We are not aware of any 
procedure that ensures that the estimator of a covariance matrix is also positive- 
definite. To overcome this difficulty, we considered the observable O = {£ — ps) 2 , 

7 We could use more general cutoff functions, but this rectangular cutoff is the most convenient 
for the present purposes. 



This has been found empirically to be true in the SW algorithm for the 2D Ising 



Potts models JL3|, and is confirmed here for the 3-state Potts model (see Section 5.5 ) 



531 and 4-state 



The bias on the estimator r- m t,o also induces a bias on the estimated variance (4.3) of the 



sample mean O. This bias is of order 1/n 2 , i.e. a factor 1/n down from the variance (|4.3j) itself. 
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which can be studied using the standard method. As we do not know exactly the 
value of /if, we can use instead the sample mean £ (which should be computed first). 
To leading order in 1/n this procedure gives the right error bar for the specific heat. 

2. We have a similar problem when computing the error bar for the second- 
moment correlation length £, which is a function of the two quantities x an d F [cf. 
( p.6|) 1- In this case we considered the random variable 

M 2 T 

& = —-—, (4.16) 

which has automatically a zero mean. The error bar for the second-moment correla- 
tion length can be written easily as a function of the susceptibility Xi the quantity 
F, and the variance of the above-mentioned observable O' . With this trick we take 
into account the cross-correlation between Ai 2 and T . This method needs the mean 
values ix m 2 an d [Lt\ in practical situations they are substituted by the corresponding 
sample means Ai 2 and T (which must be computed first). 

This is the standard procedure we have used to analyze each of our Monte Carlo 
runs. Estimates coming from multiple independent runs (for L > 256) were merged 
using the standard formulae for statistically independent data. The results are re- 
ported in Table [l] (static quantities) and Table |2] (dynamic quantities). The results 
on 71nt,£ for 16 < L < 256 can be compared with those reported in IfTIJI ; the agreement 
is excellent (x 2 = 6.17, 5 degrees of freedom, confidence level = 56%). 



4.2 Method based on independent bunches 

Unfortunately, this standard method does not provide an easy way to compute 
the error bar for complicated quantities such as the ratio t^s/Ch, which will play a 
central role in our analysis of the sharpness of the Li-Sokal bound (see Section |5.3| ) . 
We can, of course, give an upper bound on the actual error bar by using the triangle 
inequality; but this upper bound will be a significant overestimate of the true value. 
If we were then to use this overestimate in fits, we would find artificially small values 
of x 2 ; as a result, the confidence levels would be artificially high, and useless for 
distinguishing good from bad fits. (At best we could distinguish better versus worse 
fits, by looking at the relative values of % 2 .) 

This fact motivated us to look for an alternative method to compute the error bars. 
There is also another advantage in having an alternate method: we can independently 
check the assumptions and approximations made in the standard procedure. 

The second method works as follows. First, we split the whole sample of n MC 
measurements {0%, 2 , ■ ■ ■ , O n } into m bunches of equal (or almost equal) length 
£ = n/m. For each bunch % (1 < % < m) we can compute the sample mean of the 
observable O and an estimate of its integrated autocorrelation time (Oi and Ti n t,o,i-, 
respectively). Indeed, we can also estimate the corresponding variances, but these 
estimates do not play any role in what follows. When computing the autocorrelation 
functions within each bunch, we used the whole sample mean O as our estimate of fio 
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(instead of the bunch sample mean 0j); this trick reduces the bias in the estimates 
of the autocorrelation functions. 

In this way we obtain two sequences of single-bunch estimates [Oi] and {Tint,©,*}- 
If the bunches are long enough (i.e. £ ^> r int ), then the estimates from distinct 
bunches are almost statistically independent. Thus, we can define our estimates as 
follows: 

1 m 

& = -Y^Oi (4.17a) 



m *=i 



in 



vaf(a') = 1 - YsiOt-O'Y (4.17b) 
mm — 1) . 



i=i 



'int.O 



v ar(r into J 



-1 m 

(4.17c) 



m *=i 



mm 



^E^nt^-^t.o) 2 (4-17d) 



The quality of the results depends on the total number of measurements n and the 
number m of bunches we use. The merging of data coming from different runs is 
trivial in this case. 

For 4 < L < 256 we have extremely good statistics (n ~ 10 6 Ti nt) £). This allows 
us to vary m over at least one order of magnitude, and thereby to provide a cross- 
check on the standard time-series analysis. In this discussion it is useful to divide the 
observables into three categories: linear static, non-linear static, and dynamic. The 
first category includes £, £', Af, M 2 and J 7 , whose sample mean values are linear in 
the raw MC data. The second category includes the specific heat and the second- 
moment correlation length, whose mean values are non-linear functions of the raw 
MC data. Finally, the third category contains all the autocorrelation times, which 
are also non-linear functions of the raw data. 

For 4 < L < 256 we first divided the whole sample into m = 100 bunches, each 
of them with a length i ~ 10 r^s- F° r 4 < L < 64 we repeated the analysis using 
m = 1000 bunches, each of them with a length i ~ 10 3 r- mt ^. 

For the linear static observables, the "standard" and "bunch" methods always 
give identical mean values; this is a trivial identity, provided that the bunch lengths 
are exactly equal. As for the error bars on these observables, we find unsystematic 
discrepancies between the estimates given by the two methods: for m = 100 the 
discrepancies are of order 10%, and for m = 1000 they are of order 2%. In other 
words, the size of these discrepancies is roughly of order ~ 1/ \fm, and the sign is 
random; this is exactly what one expects on theoretical grounds for the statistical 
fluctuations in the estimators (4.17b) and (|4.17d| ). 



For the non-linear static observables, the agreement between the mean values 
coming from the standard and the bunch methods is exact only in the case of the 
specific heat; this is because we have used the same estimator E in both methods. For 
the correlation length, the mean values show small systematic discrepancies between 
the two methods, of order 0.05-0.1 standard deviations when m = 100; the bunch 
method always gives a larger estimate than the standard method. In absolute value, 
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these discrepancies range from 6 x 10~ 4 (for L = 4) to 2 x 10~ 2 (for L = 256). 
The same qualitative behavior is found when we repeat the analysis with m = 1000 
bunches. Now the differences in the mean value of the correlation length are 0.6-1.6 
standard deviations. In absolute value, they range from 8 x 10~ 3 (for L — 4) to 
4 x 10~ 2 (for L = 64), i.e. they are roughly one order of magnitude larger than in the 
m = 100 case. Thus, the discrepancies in the mean value of the correlation length 
are systematic with a size of order ~ m; this is exactly what one expects for the 
mean value of a non-linear observable, which is afflicted by a bias of order 1/n in 
the standard method versus l/£ = m/n in the bunch method. Regarding the error 
bars, we find unsystematic discrepancies of order 10% when m = 100, and of order 
2% when m = 1000. They behave as the error bars of the static linear observables. 

For the autocorrelation-time estimators, we find systematic differences between 
the two methods: the estimates coming from the standard method are consistently 
smaller than those coming from the bunch method. For m = 100, these differences 
are rather small compared to the statistical error bars (^ 0.5 standard deviations). 
When m = 1000 these discrepancies are much more relevant: their size is roughly 
two standard deviations. In absolute value, the discrepancies when m = 1000 are one 
order of magnitude larger than when m = 100. Thus, we find a systematic bias of size 
~ m, again as expected for a quantity which exhibits a bias of order 1/n versus l/£. 
The error bars for the autocorrelation times are consistently larger for the standard 
method than for the bunch method, except for the specific heat where the behavior is 
consistently the opposite one. For m = 100 bunches, these discrepancies are of order 
15%; and for m = 1000 the size remains at the same level. Thus, the discrepancies 
for the error bars of dynamic observables do not depend much on m; rather they are 
of order ~ 1. 

From the above discussion, we conclude that for the linear static observables the 
two methods show excellent agreement, both for the mean values (trivial equality) and 
for the error bars (random discrepancies of order 1 / \fm) . The same holds for the error 
bars of the non-linear static observables. This confirms that the standard method is 
giving accurate estimates of the error bars, at least when the run length n is large 
enough to provide a good determination of the autocorrelation time (which largely 
determines the static-quantity error bar). This l/\fm dependence also confirms our 
theoretical prediction that the standard method gives a more accurate estimate of 
the error bars than the bunch method. Indeed, the bunch method can be considered 
roughly equivalent to employing the standard method with the window width M taken 
of order the bunch length l\ but if I ^> v m t,o (as it must be in the bunch method), 
this is an unnecessarily large window width, and thus leads to unnecessarily large 
statistical fluctuations. 

For the mean values of the non-linear observables (both static and dynamic), 
we likewise confirm the validity of the standard method. Once again the standard 
method is more reliable than the bunch method, but for a different reason: the bias of 
order 1/n is much smaller than the bias of order 1/ 1 — m/n. The latter bias becomes 
particularly serious when the number m of bunches is large (as it must be in order 
to get good estimates of the error barsl). 

Finally, the least understood piece is the determination of the error bar of the 
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autocorrelation time: our data from the bunch method suggest that the standard 
method may be making a systematic error of order 15%. Perhaps this systematic error 
(if indeed it is real) arises from our neglect of the contributions of the fourth-order 
cumulant k to the covariance ( |4.10|) . This point definitely merits further investigation. 

In summary, we think that the bunch method provides a good confirmation of 
the estimates given by the standard method. We shall hereafter consider the values 
given by the standard method to be the definitive ones, except for the ratio r inti£ / Ch 
where the standard method does not yield any correct error bar. In this latter case 
we shall use the central value coming from the standard method (which in fact agrees 
with the bunch-method value to within 0.1-1%; the discrepancies are of order 0.2 
standard deviations), and the error bars coming from the bunch method with 100 
bunches (for 4 < L < 256), 26 bunches (for L = 512) or 55 bunches (for L = 1024). 
We shall also compute the upper bound on the error bar coming from the standard 
method combined with the triangle inequality. These results for r int £ /C^ are shown 
in Table ^. 



5 Data analysis 

For each quantity O, we carry out a fit to the power-law Ansatz O = ALP using 
the standard weighted least-squares method. As a precaution against corrections to 
scaling, we impose a lower cutoff L > L min on the data points admitted in the fit, 
and we study systematically the effects of varying L min on the estimates A and p and 
on the x 2 value. In general, our preferred fit corresponds to the smallest L min for 
which the goodness of fit is reasonable (e.g., the confidence levelQ is > 10-20%), and 
for which subsequent increases in L min do not cause the % 2 to drop vastly more than 
one unit per degree of freedom. 

Our final estimates for static and dynamic critical exponents are collected in Table 



5.1 Static quantities 

There are a few exactly known results concerning the 2D 3-state Potts model. We 
know all the critical exponents |TP|, |2T| and, in particular, the ratios 



1 = — » 1.73333 (5.1) 
v 15 v 1 



10 "Confidence level" is the probability that \ 2 would exceed the observed value, assuming that 
the underlying statistical model is correct. An unusually low confidence level (e.g., less than 5%) 
thus suggests that the underlying statistical model is incorrect — the most likely cause of which 
would be corrections to scaling. 



18 



which are the relevant quantities we can directly estimate from our Monte Carlo data. 
The leading correction-to-scaling exponent is also known [34, |35|: 



e = ~ , (5.3) 

using the notation in which the correction term is \j3 — (3 c \ e or L~°/ u . 

We can write the singular part of the free energy as a function of the thermal field 
t — (3 — 0c, the ordering field h, the leading irrelevant field u, and the linear size of 



the system L as 36 



f a (t,h,L) = L~ d F (tL 1/u ,hL d - p/ \uL- 9/u ) , (5.4) 



where d is the dimensionality of the lattice. If we differentiate ( |5.4p twice with respect 
to the thermal field t and then take the limit t = h = 0, we get the specific heat at 
criticality on a finite lattice: 

C H (0,0,L) ~ -^(0,0, L) = L a ' v G (0,0, uL' 6 ^) (5.5a) 

= L a/v [A + BL- e/v + ■ • •] , (5.5b) 

where a/v = 2/u—d, G is the second derivative of F with respect to its first argument, 
and the dots indicate subdominant corrections. Thus, the corrections to the specific 
heat at criticality are given by L~ A with A = 9/u = 4/5 = 0.8. A similar analysis can 
be carried out for the magnetization and the susceptibility, giving again corrections 
proportional to L~ A . 

The energy E is obtained by differentiating of the full free energy f — f s + f ns with 
respect to the thermal field t. The contribution of the non-singular piece is believed to 
be trivial: there is numerical evidence that f ns (t, L) = f ns (t, oo) |3"6| . In other words, 
this contribution has no L-dependence, and gives merely the infinite-volume value 
of the energy at the given temperature, E((3, oo). The contribution of the singular 
piece can be obtained by differentiating ( |5.4j) once with respect to t; it goes to zero 
like L~ d+1 / U , which thus gives the leading correction to scaling for the energy. This 
correction is of order L -4 / 5 , which (by pure coincidence as far as we can tell) is exactly 
the same order as the correction L~ A for the divergent static observables. Finally, we 
note that the energy of the 2D q = 3 Potts model at criticality is also exactly known 
II to be E(p c , oo) = 1 + l/y/3 w 1.577350. 

We can check these predictions by performing the fit E — E(/3 C , oo) = AL~ W . The 
quality of the fit is very good already for L min = 16: 

w = 0.803 ± 0.002 (5.6) 

with x 2 — 3.23 (5 DF, level = 66%). The agreement with the predicted exponent 
2 — 1/v = 4/5 = 0.8 is truly spectacular. (Indeed, it is probably a coincidence that 
the agreement is so good. In general, one can't expect to obtain anywhere near this 
accuracy for correction-to-scaling exponents.) 
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The fits of the susceptibility to a power law AL^ U are quite stable. Our preferred 
fit corresponds to L min = 64: 

^ = 1.73444 ± 0.00043 (5.7) 

with x 2 = 4.09 (3 DF, confidence level = 25%). This result is 2.6 standard deviations 
away from the exact value 7/z/ — 1.73333. This discrepancy could be due to correc- 
tions to scaling. We can try to fit our data to AL 26 / 15 (1 + BL~ A ) with various choices 
for A = 1.1, 1.0, . . . ,0.1 as well as x log (i.e., a correction 1/logL). We find that 
the best fits correspond to A « 0.8, in agreement with the theoretical prediction; for 
A = 0.8 and L min = 8 we obtain x 2 — 6.55 (6 DF, level = 36%). Surprisingly, in all 
these fits we find that the x 2 remains almost constant when L min is increased beyond 
our "preferred" value (and the confidence level consequently deteriorates): for exam- 
ple, with L min = 256 we get x 2 — 3.37 (1 DF, level = 7%) for the pure power-law 
behavior and x 2 = 3.91 (1 DF, level = 5%) for the fit with 7/z/ = 26/15 and A = 0.8. 
This suggests that the point with L = 1024 is off by about two standard deviations 
(possibly because the error bar is underestimated^). As a matter of fact, if we drop 
this point, we obtain a good power-law fit for L min = 128: 

- = 1.73337 ±0.00080 (5.8) 

with x 2 — 0-39 (1 DF, level = 53%). This result agrees excellently with the exact 
value. If we impose the right power 7/z/ = 26/15 and try to fit the first correction- 
to-scaling exponent A, we again find that the best fits correspond to A rs 0.8. For 
Lmin = 8 we get x 2 = 3.84 (5 DF, level = 57%). 

For the specific heat we find that the fits to power law AL a l v are not stable at all: 
the confidence levels are horrible, and there is a clear trend towards smaller values of 
ajv as L min is increased. The least bad fit is obtained for L min = 256: 

- = 0.4240 ± 0.0030 (5.9) 
v 

with x 2 — 3-80 (1 DF, level = 5%). This value is eight standard deviations away 
from the exact result ajv = 2/5 = 0.4. Unlike the 4-state Potts model 24|, we 
do not expect multiplicative corrections to the leading term of the specific heat. We 
do, however, expect additive corrections to scaling of the form AL 2 ^{\ + BL~ A ). If 
we try the same exponents A as in the susceptibility, we find a decent fit for A r* 0.6: 
for L min = 128 we get x 2 — 1-21 (2 DF, level = 55%). This value of the exponent A 
is not far from the expected value 0.8, but the cause of the discrepancy is unknown; 
perhaps there is a large next-to-leading correction to scaling. 

Finally, the second-moment correlation length £ is expected to behave linearly 
in L as L — > 00. In particular, the ratio x = £/L should approach a constant x*. 

11 This is quite possible: though the total run length at L = 1024 is 7 x 10 4 Ti n t,£, the individual 
runs (on which the time-series analysis was performed) ranged in length from only 10 Ti n t,£ to 
2.5 x 10 Tint g; and with such short runs the time-scries analysis may not be completely reliable. 
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We have tested this behavior. Already for L m i n = 16 our data are consistent with a 
constant value 

x* = lim = 0.93235 ± 0.00033 (5.10) 

with x 2 — 8-08 (6 DF, level = 23%). The fact that a good fit can be obtained 
with such a small L min implies that the corrections to scaling are very small for this 
observable. [If we fit to £{L)/L = x* + BL~ A with 0.1 < A < 1.5, we find a very slight 
improvement in the goodness of fit, and the estimated x* decreases somewhat (by less 
than 0.0008 if A > 0.6). But the estimated coefficient B is consistent with zero at 
the 1.5er level.] As in the case of the susceptibility, we also find that the goodness of 
fit deteriorates as L min is increased (e.g., for L min = 512 we have \ 2 — 4.80, 1 DF, 
level = 3%) This might be due to the fact that the value for L = 1024 is a little bit 
off (or its error bar is underestimated). If we drop this point, we obtain a good fit 
again for the same L min = 16: 

x* = 0.93229 ±0.00034 (5.11) 

with x 2 — 5.82 (5 DF, level = 32%). But now the fit with L m i n = 256 is reasonable 
(x 2 = 1-24, 1 DF, level = 27%). We remark that the value x* should in principle 
be calculable by conformal field theory; we hope that someone will perform this 
calculation and test our results ( |5.10| )/( |5.11| ). 



5.2 Dynamic quantities 

In this section we are going to fit the autocorrelation times for the observables 
O = S, £', N and M 2 to a simple power law r intiC > = AL Zint <° . 

Let us start with the energy £. The fit T- mt) £ = AL Zint ' £ is not very stable: the 
estimate of the power decreases systematically as L min is increased, and the x 2 is 
poor until L m i n = 128 where the estimate stabilizes within errors (this behavior 
suggests that there are strong corrections to scaling). Our preferred fit corresponds 
to L min = 128: 

Zint,s = 0.515 ±0.006 (5.12) 

with x 2 — 0.44 (2 DF, level = 80%). This value is greater than our estimate for a/u; 
thus, the Li-Sokal bound ( |1.2|) holds, though apparently not as a strict equality. To 
compare our result with the one reported in [TT], we redid the fit using only the data 



with L < 256. The fit again shows a systematic decrease of the exponent as L min is 
increased, as well as a poor x 2 , until L min = 64 where we get z intt £ = 0.531±0.005 with 
X 2 = 2.67 (1 DF, level = 10%). This is consistent with the result z int:£ = 0.55 ± 0.02 
reported in |[Llj| , but our error bar is one-fourth of theirs. The slightly higher estimate 
of Zint,£ in |TT[] seems to arise from corrections to scaling induced by their choice 



Lmin = 16. Actually, the best fit to the data in |TT| corresponds to L min = 32 and 
gives z int>£ = 0.52 ± 0.02 with X 2 = 3.93 (2 DF, level = 14%). 

The fit of the autocorrelation time for the bond occupation M follows the same 
pattern: the estimate of the power decreases strongly as L min increases, and the x 2 
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is initially horrible; eventually the power stabilizes within errors, and the x 2 becomes 
reasonable. Our preferred fit is L min = 128: 

z in t,Af = 0.529 ±0.006 (5.13) 

with x 2 — 0.71 (2 DF, level = 70%). The difference with respect to z intt s is only 2.3 
standard deviations, consistent with the theoretical prediction (|2.29|) that z int ^ = 



z int,e- To check this result, we studied the ratio t^^/t^s- Since the standard time- 
series-analysis method would give only a upper bound on the right error bar for this 
quantity, we used instead the error bar provided by the "bunch" method. Among 
all the Ansatze we tried, only one gave a good fit for L min = 128: asymptotically 
constant with corrections ~ L -1 / 4 (x 2 = 0.75, 2 DF, level = 69%). Thus, we conclude 
that in the SW algorithm the dynamic critical behavior of the energy and the bond 
density are the same: z int ^ = z int ^. 

In the same way, we considered the nearest-neighbor connectivity £'. The pattern 
of the fit is the same as above, and our preferred fit again corresponds to L min = 128: 

z in t,e> = 0.514 ±0.006 (5.14) 

with x 2 — 0.47 (2 DF, level = 79%). This time, the agreement with z- mtj s is extremely 
good, confirming the theoretical prediction ( |2.41| ) that z int)£ > = z intt £. 



Finally, a similar behavior is observed in the fit for the autocorrelation time of the 
squared magnetization M 2 : the estimate of the power shows a clear trend towards 
smaller values as L min is increased, and the x 2 is initially poor. Our preferred fit is 
Lmin = 128, for which we get 

z . mt M2 = 0.475 ±0.006 (5.15) 

with x 2 — 0.36 (2 DF, level = 84%). This power is slightly smaller than z- mt £ (the 
difference is seven standard deviations). It would be interesting to know whether this 
difference is significant or not. Let us consider the ratio r int t M 2 / T mt,s, again using 
the error bar provided by the "bunch" method. We tried to fit this ratio to various 
Ansatze, but only two gave good results for L min = 32: a pure power-law behavior 
ALP with p = -0.0409 ± 0.0007 (x 2 = 0.57, 5 DF, level = 97%); and a constant plus 
corrections of the type ~ L" 1 / 16 (x 2 = 0.41, 5 DF, level = 98%). We are therefore 
unable to resolve whether z- mtM 2 is exactly equal to z int ^ or not; but if it is not equal, 
then it is only very slightly smaller (z- mtjM 2 — z int)£ w —0.04). 

In conclusion, we have shown numerically that the energy S, the bond occupation 
A/", and the nearest-neighbor connectivity £' all have the same dynamic critical ex- 
ponent 2i n t, as was "almost proved" in Section ^. On the other hand, the observable 
M 2 has a similar but perhaps not identical dynamic critical behavior: z^^ may 
coincide with z int ^, or it may be slightly smaller. 

5.3 Analysis of the Li— Sokal bound 

In the previous subsection we have seen that all the observables considered have 
the same (or, in the case of A^ 2 , almost the same) dynamic critical exponent, and the 
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common z- mt is strictly larger than a/u. This implies that the Li-Sokal bound ( |1.2j ) 
is not sharp. However, the difference z int — a/u ~ 0.115 is actually not very large. 
There are a few arguments in favor of a more detailed analysis: 

i) The power-law fit to the specific heat was not very good: the estimated value of 
a/u decreased as L min increased, and we were unable to reach (within statistical 
errors) the exact value a/u = 2/5. If we compare the observed values of z int 
and a/u, we find that z- mt — a/u is only w 0.09. 

ii) The power-law fits to the autocorrelation times also exhibited this trend to 
smaller values as L min increased. Even though the fits seemed stable for L min = 
128, this might well be an artifact due to the large error bars associated with 
the largest lattices {L > 512) compared to the smallest ones (L < 256). 



iii) In [12, [13] it was shown that differences z int — a/u of order 0.1 could actually 



be due to multiplicative logarithmic corrections. Indeed, such a multiplicative 
logarithmic correction was found to be a likely scenario even for models not 
having such a logarithmic correction in the specific heat. 



Here, we will follow the approach of Refs. |L2], |13[, and consider the ratio v m t,e/CH\ 
we shall use the error bars obtained from the "bunch" method (see Table 0)-0 We 
have tried to fit this ratio to different Ansatze: 



1) A pure power- law behavior AL P . 

2) A logarithmic growth, either as A + BlogL or as A\og p L. 

3) Asymptotically constant with additive corrections to scaling A + BL~ A . We 
have considered the cases A = 2, 1, ~, |, |, and x log (i.e., t-^^/Ch = A + 
B/logL). 

Among all these Ansatze, only two were reasonably good. The best one is the 
simple power-law behavior ALP . For L m i n = 32 we obtain 

p = 0.084 ±0.002 (5.16) 

with x 2 = 1-72 (4 DF, level = 79%). When L min > 32, the value of the power p stays 
stable, and the x 2 decreases slowly and consistently. 

The second-best fit corresponds to the logarithmic growth A + B log L. For L min = 
64, we get x 2 — 1-93 (3 DF, level = 59%). Again, the estimates are stable for 
Lmin > 64 and the x 2 is reasonable. However, this Ansatz seems to be slightly 
inferior to the power-law fit: both L m i n and the x 2 are greater than in the power- law 
case. To test the logarithmic Ansatz, we can impose the known critical behavior of 
the specific heat and perform the fit T- mtj£ = L 2 ^ 5 (A + BlogL). A reasonably good 
result is obtained for L min = 32, giving x 2 — 1-54 (4 DF, level = 82%). 

12 In this section we only consider the energy E, as the other observables have the same (or very 
slightly smaller) dynamic critical exponent. 
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The asymptotically constant fits are always horrible, unless one take L m i n — 256 
or larger. The only semi-exception is A = |, which gives a tolerable fit (x 2 = 5.54, 3 
DF, level = 14%) already for L min = 64. 

Of course for L min = 256 we obtain reasonably good fits for all Ansatze. But this 
is because the error bars for L > 256 are so large that we are unable to distinguish 
between the various scenarios. 

In summary, we have shown that there are only two likely scenarios for the ratio 
r mt,s/ Ch- 

\ AL P with p = 0.084 ± 0.002 , r _ 
n«,e/C H ~{ a + B]orL P (5-17) 



or equivalently 



f AL Z ^> S with z int g = 0.515 ± 0.006 , , 

~ \ L^{A + B\ogL) (5 ' 18) 



The first scenario implies that the Li-Sokal bound (|1.2|) is a strict inequality, with 
p = Zinb t g — ajv « 0.08-0.12; while the second scenario means that this bound is 
violated only by a logarithm (i.e., the bound is sharp modulo logarithms). Thus, the 
same "dynamic continuity" we found interpolating between the 2-state and the 4-state 
Potts models along the AT self-dual line is also found along the g-state Potts model 
critical line. Moreover, the power p found is midway between the power found for the 
Ising model (p ~ 0.05) and the one found for the 4-state Potts model (p ~ 0.12) [|13j. 

5.4 Further discussion of the Li— Sokal bound 



The proof of the Li-Sokal bound (Section |2.2|) is based on computing the au- 
tocorrelation function for the observable Af at time lags and 1 in terms of static 
observables, and then exploiting the inequality pjvov(£) > Paov(1)'*'- The apparent 
non-sharpness of this bound indicates that the large-t behavior of PMj^it) is not fully 
predicted by its behavior at t = 1. Otherwise put, if we define the initial autocorre- 
lation time 

1 l+poo(l) /, 10 s 

2 l-poo(l) 

then the Li-Sokal bound 

Thbtf > 7init,A/" = z ft + o (5.20) 

1 — p hi L 

is apparently not sharp: that is, r int ^ and r init ^ diverge with different critical expo- 
nents z intj M and ZMt^r = & jv. 

One might ask whether the situation could be improved by computing Cj\fx(t) ex- 
actly at (for example) time lag t = 2 or t = 3. Because of the identities fl2.25|) /( [2.36|) , 
this would be equivalent to carrying out the Li-Sokal proof using S or £' as the test 
observable in place of A/". At present we do not know how carry out this computa- 
tion — the trouble is that we do not know how to express (7^7^) in terms of spin 
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observables — but we can nevertheless test numerically whether we would obtain a 
better bound on Zi n t,e and/or z- int>e i if we could carry out this computation. 

We thus computed r initjC) for O = Af, £, £', with error bars given by (flKjQ we 
then tried various fits, among others the fit T init0 /CH = AL r .f^\ 

For O = M we have the identity 



r init,./V 



p 



1 



1 - p E ' 2C H 
p 1 

1 - pE((3, oo) 



+ aL~^ + 



bL- 2 / 5 + cL~^ + ■ 



(5.21a) 
,(5.21b) 



where E((3,oo) is the infinite-volume energy at inverse temperature (3; and we have 
taken into account the leading terms and the corrections to scaling discussed in Sec- 
tion |5TI| . At criticality (3 = (3 C = log(l + a/3) [hence p = p c = a/3/ (1 + \/3)] we know 
fill the exact constant term: 



Pc 



1 - p c E(/3 C , oo) 1 + a/3 



1.098076 . 



(5.22) 



Indeed, if we fit our data to the Ansatz t^^m/Ch — 3/(1 + a/3) — bL s we obtain a 
good fit for L min = 16, and the estimate is 



0.415 ±0.017 



(5.23) 



(x 2 = 1-92, 3 DF, level = 86%), in excellent agreement with the theoretical prediction. 
Of course, we could try also the more primitive fit T init A r/C^ = AL r , and we get a 
good fit for L min = 128: 

r = -0.010 ±0.002 (5.24) 

(x 2 = 0.90, 2 DF, level = 64%), which is very close to the exact result r = 0. All 
this is trivially to be expected, since our raw data for pj\fj\f(l), Ch and E are in good 
agreement with the theoretical identity ( ^.16| ). 

The analysis becomes non-trivial when we look at O = £ and £' . For O = £ we 
fit Tinit^/Cjy = AZ7, and get a good fit for L min = 64: 



(x 2 



0.91, 3 DF, level 



r = 0.003 ±0.003 (5.25) 
For O = £' the analogous fit is good for L min = 32: 
r = 0.016 ±0.001 (5.26) 



(x 2 — 0.51, 4 DF, level = 97%). These exponents r are a factor of ~ 5 smaller 
than the exponent estimates for p = z — a/u found in Section [5~3| . This suggests 

13 Since the autocorrelation fun ction for these observables is almost exactly a pure exponential, 
it suffices to perform the sums in ( 4.10 ) analytically and then insert the appropriate value of r (see 
11). 

1 For the error bars on Tini^o/C/f, we used for simplicity the triangle inequality rather than the 
(more correct) bunch method. 
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that Ti n it,£ / Ch and T in i ti £' / Ch in fact tend to finite constants as L — > oo (as we know 
rigorously to be the case for t^^/Ch), and that the apparent non-sharpness of the 
Li-Sokal bound cannot be remedied by using £ or £' as a test observable in place of 
N '. Rather, the non-sharpness of the bound arises from the fact that the long-time 
behavior of pj^j^(t) is not sufficiently well predicted by its behavior at any small time. 

5.5 Exponential autocorrelation time 

The exponential autocorrelation time for an observable O is defined as0 

Te W ,o = lim ^L- . (5.27) 

t^oo log p OQ (t) 

This autocorrelation time measures the decay rate of the "slowest mode" of the sys- 
tem, provided that this mode is not orthogonal to O. 

The critical behavior of r exP) o is, in general, different from the behavior of T intj e>. 
This fact can be seen from the standard dynamic finite-size-scaling Ansatz for the 
autocorrelation function poo{t): 

P oo(t;L) » \ t \-voho(—;^) ■ (5.28) 

(Here the dependence on the coupling constants has been suppressed for notational 
simplicity.) Summing (|5.28j) over t, it follows that 

r^o ~ r^ P S , (5.29) 

or equivalently, 

Zmt,0 = (1 -Po)Zexp,0 ■ (5-30) 

Thus, only when p = do we have z int: o = ^exp,o @ • I n this latter case the Ansatz 
( 5.28Q can be rewritten in the equivalent form 



P oo(t;L) « U—;^-) • (5.31) 

To test this latter Ansatz, we have plotted logpoo(t) versus tjr^o for the ob- 
servables O = Af (Figure [l]), S (Figure §) and £' (Figure §). On each figure we have 
plotted the data coming from different lattice sizes (4 < L < 128) with different 
symbols; the error bars are computed from (|4.10|) , using for simplicity the approxi- 
mation [|l^] that the decay is a pure exponential (which is here almost exact). On 
each graph, we have also depicted for reference a line corresponding to the pure ex- 
ponential poo{t) = e - '/ Tint <°. In these plots we have omitted the data for L > 256 for 

15 For a general Markov chain, the "lim" should strictly speaking be replaced by "lim sup", and 
Poo(t) should be replaced by its absolute value. But here it can be proven that the limit really 



exists, and that poo(t) > for all t; this follows from the spectral representation ( 2.19 ) [or rather 
its analogue for O] . 
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the sake of visual clarity; these data agree well with the curve found for smaller L, 
but with huge statistical error bars. 

For O = M we see that the data coming from 16 < L < 128 collapse well onto a 
single curve. The lattices L = 4, 8 show slight systematic deviations from this limiting 
curve: these deviations are negative for t/T int ^ < 1.5 and positive for t/r int ^ > 1.5. 
This trend continues for 16 < L < 128, but the deviations are in most cases smaller 
than the error bars (especially for the larger lattices). 

A similar behavior is found for O = £ (Figure 0) and O — £' (Figure but the 
corrections to scaling are much weaker, and their sign is opposite from those seen for 
M. 

Thus, within statistical errors, we have found that the Ansatz ( |5.31|) is satisfied. 
This implies that the integrated and exponential autocorrelation times for these three 
observables have exactly the same dynamic critical exponent, i.e. z- m tp — z e*p,o f° r 
O = Af, £, £'. This equality does not hold as a general rule in the theory of dynamic 
critical phenomena 0, ||, but it does appear to hold for algorithms of Swendsen-Wang 
type. 

Finally, we find that pj^_\r(t) differs slightly but noticeably from the pure expo- 
nential e - */ 7 " 111 '-^. The discrepancy from a pure exponential becomes smaller when we 
consider the other two observables £ and £'. This is to be expected theoretically from 
the relations £ ~ PbondA/" and £' ~ -P sp in£ [cf. (|2.14j )/ (|2.35|) 1: each action of Pbond or 



P sp in helps to "purify" the slowest mode, so that the autocorrelation function becomes 
closer to a pure exponential. On the other hand, the identities ( |2.26|) / Q2.37|) imply^] 
that the limiting (scaling) functions h_\f, and are identical (assuming they exist 
at all) ; this means that at least some of the curves in Figures [l]-Q have not yet reached 
their scaling limit. 

Finally, a crude fit suggests that r- int ^/r cxv , y s w 0.96, in agreement with the idea 
that pes(t) is almost but not quite a pure exponential. 
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Using the fact that puuiX) —* 1 and Tmt,AA — ► oo as L — > oo, as follows from the Li-Sokal 



identity ( 2.16 ) and the fact that the specific heat Ch is divergent. 
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L 


MCS 


X 






E 


4 


exact 


12.204711 


1.496719 


3.862380 


1.710062 


4 
8 
16 
32 
64 
128 
256 
512 
1024 


4.9 
6.9 
9 9 
14.9 
19.9 
29.9 
40.8 
12.9 
5.5 


12.2093 ± 0.0054 
41.3527 ± 0.0204 
138 4491 ± 0723 
462.7014 ± 0.2408 
1542.6921 ± 0.8428 
5135.9512 ± 2.7481 
17082.8221 ± 9.2802 
56760.2838 ± 65.0770 
189676.5530 ±387.3500 


1.4952 ±0.0021 
2.4712 ±0.0031 
3 71 62 + 0046 
5.2941 ±0.0064 
7.3839 ±0.0094 
10.1007±0.0127 
13.7116 ±0.0176 
18.4752 ±0.0507 
24.5281 ±0.1247 


3.8635 ±0.0043 
7.5419 ±0.0072 
14 9267 ±0 0134 
29.8509 ±0.0250 
59.6964 ±0.0501 
119.3401 ±0.0952 
238.4565 ±0.1888 
475.9494 ±0.7779 
958.8929 ±2.7875 


1.71051 ±0.00039 
1.65391 ±0.00026 
1 62070 ±0 00016 
1.60223 ±0.00010 
1.59163 ±0.00006 
1.58552 ±0.00003 
1.58202 ±0.00002 
1.58001 ±0.00003 
1.57893 ±0.00003 


oo 


exact 








1.577350 



Table 1: Static data from the MC simulations at the critical point of the 3-state Potts 
model. For each lattice size (L), we include the number of performed measurements 
(MCS) in units of 10 6 , the susceptibility (x), the specific heat (Ch), the second- 
moment correlation length (£), and the energy (E). The quoted errors correspond to 
one standard deviation (i.e. confidence level ~ 68%). The first row ("exact") gives 
the exact results for the 4x4 lattice, and the last row ("L = oo") gives the exact 
energy in the limit L — > oo. 



L 


T int,M 2 


r mt,Af 


r int,£: 




4 


4.013 ±0.018 


3.205 ±0.013 


4.023 ±0.018 


4.034 ±0.018 


8 


5.983 ±0.028 


5.194 ±0.023 


6.033 ±0.028 


6.080 ±0.028 


16 


8.816 ±0.041 


8.084 ±0.036 


9.025 ±0.043 


9.117±0.043 


32 


12.648 ±0.057 


12.200 ±0.055 


13.280 ±0.062 


13.438 ±0.063 


64 


18.101 ±0.085 


18.313 ±0.086 


19.549 ±0.095 


19.769 ±0.097 


128 


25.665±0.117 


27.094 ±0.127 


28.525 ±0.137 


28.820 ±0.139 


256 


35.691 ±0.164 


39.204 ±0.189 


40.824 ±0.200 


41.212 ±0.203 


512 


49.775 ±0.480 


56.667 ±0.583 


58.511 ±0.612 


58.974 ±0.619 


1024 


68.387± 1.183 


80.489 ±1.511 


82.488 ±1.567 


83.030 ±1.583 



Table 2: Autocorrelation times for the runs performed at the critical point of the 3- 
state Potts model. For each lattice size (L), we include the integrated autocorrelation 
times for the squared magnetization (r int ^2), the bond occupation (r inti Ar), the energy 
( r mt,£), and the nearest-neighbor connectivity (r inti £'). 
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L 


r int,,f/ Ch 


4 


2.6906 


± 


(0.0111, 


< 


0.0161) 


o 




2.4413 


± 


(0.0092, 


< 


0.0145) 


16 


2.4285 


± 


(0.0090, 


< 


0.0145) 


32 


2.5084 


± 


(0.0084, 


< 


0.0147) 


64 


2.6476 


± 


(0.0111, 


< 


0.0163) 


128 


2.8241 


± 


(0.0111, 


< 


0.0171) 


256 


2.9773 


± 


(0.0113, 


< 


0.0184) 


512 


3.1670 


± 


(0.0281, 


< 


0.0418) 


1024 


3.3630 


± 


(0.0442, 


< 


0.0810) 



Table 3: Ratios t^^/Ch for the runs performed at the critical point of the 3-state 
Potts model. In parentheses we give our error-bar estimates. The first number shows 
the error bar coming from performing the "bunch" method with 100 bunches. The 
second number is obtained by using the triangular inequality with the numerical 
results coming from the standard method (see Section £|). 



Exponent 


numerical exact 


ajv 


1.73444 ±0.00043 26/15 
0.4240 ±0.0030 2/5 


&«&,£' 
z mt,M 2 


0.515 ±0.006 > 2/5 
0.514 ±0.006 > 2/5 
0.529 ±0.006 > 2/5 
0.475 ± 0.006 



Table 4: Numerical estimates for the static and dynamic critical exponents of the 3- 
state Potts model (second column). We also include the exact results for comparison 
(third column). 
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Figure 1: Plot of pjvovW versus t/r hlt ^ for 4 < L < 128. The different symbols 
denote the different lattice sizes: L = 4 (*), L = 8 (+), L = 16 (x), L = 32 (□), 
L = 64 (O), and L = 128 (o). We have also depicted the line corresponding to the 
pure exponential puuif) = exp(-t/r int>A r). 
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Autocorrelation Function of E 




/ T int,E 



Figure 2: Plot of pee(t) versus t/r int)£ for 4 < L < 128. The symbols are as 
in Figure [I]. We have also depicted the line corresponding to the pure exponential 
Pseif) = exp(-t/r int>£ -). 
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Autocorrelation Function of E' 




Figure 3: Plot of pcs'it) versus t/r^^i for 4 < L < 128. The symbols are as 
in Figure [I]. We have also depicted the line corresponding to the pure exponential 
P£'S'{t) = exp(-t/r inti£ /). 
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